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Abstract We derive a compositional compressible two-phase, liquid and gas, flow model 
for numerical simulations of hydrogen migration in deep geological repository for radioac- 
tive waste. This model includes capillary effects and the gas high diffusivity. Moreover, it 
is written in variables (total hydrogen mass density and liquid pressure) chosen in order to 
be consistent with gas appearance or disappearance. We discuss the well possedness of this 
model and give some computational evidences of its adequacy to simulate gas generation in 
a water saturated repository. 
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1 Introduction 

The simultaneous flow of immiscible fluids in porous media occurs in a wide variety of ap- 
plications. The most concentrated research in the field of multiphase flows over the past four 
decades has focused on unsaturated groundwater flows, and flows in underground petroleum 
reservoirs. Most recently, multiphase flows have generated serious interest among engineers 
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concerned with deep geological repository for radioactive waste. There is growing aware- 
ness that the effect of hydrogen gas generation, due to anaerobic corrosion of the steel en- 
gineered barriers (carbon steel overpack and stainless steel envelope) of radioactive waste 
packages, can affect all the functions allocated to the canisters, waste forms, buffers, back- 
fill, host rock. Host rock safety function may be threaten by overpressurisation leading to 
opening fractures of the host rock, inducing groundwater flow and transport of radionuclides. 

The equations governing these flows are inherently nonlinear, and the geometries and 
material properties characterizing many problems in many applications can be quite irregu- 
lar and contrasted. As a result, numerical simulation often offers the only viable approach to 
the mathematical modeling of multiphase flows. In nuclear waste management, the migra- 
tion of gas through the near field environment and the host rock, involves two components, 
water and pure hydrogen H2; and two phases "liquid" and "gas". It is then not clear if con- 
ventional models, like for instance the Black-Oil model, used in petroleum or groundwater 
engineering are still valid for such a situation. Our ability to understand and predict under- 
ground gas migration is crucial to the design of reliable waste storages. This is a fairly new 
frontier in multiphase porous-media flows, and again the inherent complexity of the physics 
leads to governing equations for which the only practical way to produce solutions may be 
numerical simulation. This exposition provides an overview of the types of standard models 
that are used in the field of compressible multiphase multispecies flows in porous media 
and includes discussions of the problems coming from H2 gas being one of the components. 
Finally the paper also addresses one of the outstanding physical and mathematical problems 
in multiphase flow simulation: the appearance disappearance of one of the phases, leading 
to the degeneracy of the equations satisfied by the saturation. It has been seen recently in a 
benchmark organized by the French agency in charge of the Nuclear Waste management in 
France 1141 that none of the usual codes used in that field where able to simulate adequately 
the appearance or/and disappearance of one of the phases. In order to overcome this diffi- 
culty, we will discuss a formulation based on new variables which doesn't degenerate. We 
will demonstrate through two test cases, the ability of this new formulation to actually cope 
with the appearance or/and disappearance of one phase. The scope of the paper is limited to 
isothermal flows in rigid porous media and will not include the possible process of "pathway 
dilation" as described in j5J- 



2 Conceptual and mathematical model 

Our goal is first to present a survey of the conventional models used for describing two- 
phase two-components flow in porous media. Most of these models have been designed 
and widely used in petroleum engineering (see for instance: [T), (2), (3], (5), )lll). We 
consider herein a porous medium saturated with a fluid composed of 2 phases : liquid and 
gas. According to the application we have in mind, the fluid is a mixture of two components: 
water (mostly liquid) and hydrogen (H2, mostly gas). Water is present in the gas phase 
through vaporization, and hydrogen is present in the liquid phase through dissolution. The 
fluids are compressible and the model is assumed to be isothermal. For simplicity, we assume 
first the porous medium to be rigid, meaning the porosity <J> is only a function of the space 
variable <P = 4>(x); and second, we neglect the pressure induced dilation of gas pathways. 
Hydrogen being highly diffusive we have described in detail how the diffusion could be 
taken in account in the models. 
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2. 1 Petrographic and fluid properties 
2.1.1 Fluid phases 

The two phases will be denoted by indices /, for liquid and g for gas. Associated to each 
phase in the porous media are the following quantities: 

- liquid and gas phase pressures: pi, p g ; 

- liquid and gas phase saturations: Si,S g ; 

- liquid and gas phase mass densities: Pi, p g ; 

- liquid and gas phase viscosities: fJ-i, H g ; 

- phase volumetric flow rates, q/ and q g ; 

Then, the Darcy-Muskat law says: 

qi = -K(x) k -^l(V Pl -p lg ), qg = -K(x)^^(V^-p g g), (1) 
Pi Pg 

where K(x) is the absolute permeability tensor, kr\ and kr g are the relative permeability 
functions, and g is the gravity acceleration. 

The above phase mass densities and viscosities are all functions of phase pressures and 
of phase composition. Moreover, phase saturations satisfy 

S,+S g = l (2) 

and the pressures are connected through a given experimental capillary pressure law: 

Pc{S g ) = p g - Pi- (3) 

From definition ((3) we should notice that p c is strictly increasing function of gas saturation, 
p' c (S g )>0. 



2.1.2 Fluid components 

Water and pure hydrogen components will be denoted by indices w and h. Since the liquid 
phase could be composed of water and dissolved hydrogen we need to introduce the water 
mass density in the liquid phase p, w , and the hydrogen mass density in the liquid phase pf. 
Similarly, we introduce the water and hydrogen mass densities in the gas phase: p" ', p g . 
Note that the upper index is always the component index, and the lower one denotes the 
phase. We have, then 

Pi=p? + pl Pg = p g v + p' g \ (4) 

Since the composition of each phase is generally unknown we introduce the mass fraction 
of the component i e {w, h} in the phase a 6 {g, 1} 

K = ^ < + ffl£ = l,ae{g,Z}. (5) 

Pa 
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2.2 Mass conservation equation 

The conservation of the mass applies to each component i and reads, for any arbitrary control 
volume M: 

= ie{w,h}. (6) 

at 

where 

- m' is the mass of the component in the control volume 0t, at given instant t; 

- F' is the rate at which the component i is leaving (migrating from) the volume S£, at 
given instant ?; 

- is the source term (the rate at which the i component is added to M by the source). 
The mass of each component i is a sum over all phases g and /: 

m' = / <t> (S/piCo', +S g pg(o' s ,) dx, ie{w,h}. 

J & 

In each phase, the migration of a component is due to the transport by the phase velocity 
and to the molecular diffusion: 

F' = / (p/co/q, +Pg(o' g qg +jj +j g ) ■ ndx; i e {w,h}; 

where n is the unit outer normal to dS$. The phase flow velocities, q\ and q g are given by 
the Darcy-Muskat law HJ, and the component i diffusive flux in phase a is denoted j'^, 
i £ {w,h}, a 6 {g,l}, and will be defined in the next paragraph, by equations d 1 2b - From ® 
we get the differential equations: 

t>j t (S,p,w? + S g PgCOg') + div ( P/ <q, + PgCOg'qg +jf +#) = (7) 



^p/Co'^ + SgPgCo'g 1 ) +div [p^qt +p g w n g qg +j* = 9*. (8) 
2.2.1 Diffusion fluxes 

From M w and M h , the water and hydrogen molar masses, using definitions we define the 
following water and hydrogen molar concentrations in each phase a 6 {g, I}: 

,h _ SaPq _ SaPqCQq w _ S a Pq _ SqPgtoq 
Ca ~ M h M h ' C(X ~ M w M K 

Then the phase a molar concentrations, a 6 {g,l}, is: 

h , „w c „ ( w a w a \ nrt . 
c * = c a+Ca=S a pa\j^ + —). (10) 

Usually, component i diffusive flux in phase a is assumed to be depending on X' a , the com- 
ponent i molar fraction in phase a 6 {g, I}, defined from ^ and d 10b : 

j^h _ yw . V 1 X' — 1 (11) 

a ~ c a ~ co h a + (M h /M»')coS ' a ~ c a ~ coS + (M"/M>>)w* ' £ a ' 
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for a 6 {l,g}- Molar diffusive flux of component i in phase a € {g,l} is given by 

J'g = -CaD'aVX'a, a e {g,l}, i e {w,h} 

and give molar diffusive flow rate of component i through the unit area. Coefficients D 1 ^ and 
D n a (unit L 2 /T) are Darcy scale molecular diffusion coefficients of components in phase 
cc 6 {g, I}. Mass flux of component i in phase a, in equations {7), <[Sj> is then obtained by 
multiplying the above component molar diffusive fluxes, J' a , by the molar mass M; of the 
component i, and by the rock porosity: 

j h a = -&M h c a D» a VX*, r a = -4>M w c a D w a VX%. (12) 



Remark 1 : Number of unknowns in the system is eight: Si,pf ,pf,pi,S g ,p g ,p g ,p g ; but, 
up to now, we have only four equations (O, (O, ITJ and ®, and therefore, four additional 
equations are needed to close the system. 

Remark 2 : Note that D h a and D„ are not exactly molecular diffusion coefficients in phase 
a, corresponding to molecule-molecule interactions in free space but <£>D' a , a 6 {g,l}, 
i € {w, h}, are effective diffusion coefficients, obtained from the strict molecular diffusion 
coefficients by a kind of averaging through the whole porous medium (see section 2.6 in 
0, or Q3 or (6))- Moreover, for simplicity, we have not included in diffusion of compo- 
nent i in phase a any dependancy on phase saturations; and then did not consider possible 
non linear effects coming from coupling between advective-diffusive transport ("dusty gas" 
model) and molecular streaming effects (Knutsen diffusion) (see for instance 1 10]). 

Remark 3 In a binary system diffusive fluxes satisfy j* + j^j =0, for a 6 {g, /}, and there- 
fore we have 

M h D h a =M w D w a , ae{g,l}. (13) 



2.3 Phase equilibrium Black-oil model 

The additional equations needed to close the system of equations (2j, (13}, (UJ and ((8} will 
come from the assumption that the two phases are in equilibrium; equilibrium meaning 
that at any time the quantity of hydrogen dissolved in the water is maximal for the given 
pressure, and similarly, the quantity of evaporated water is maximal for the given pressure. 
Composition of each phase is then uniquely determined by its phase pressure and saturation. 
In this flow situation we say that water and gas phases are saturated. Nevertheless it could 
happen that one of the phases disappears: either water can completely evaporate or hydrogen 
can be completely dissolved in the water. Then in these situations the composition of the 
remaining phase is not uniquely determined by its phase pressure, and the phase composition 
becomes an independent variable (instead of the saturation which is now constant, or 
1). This flow situation correspond to the so-called unsaturated flow. For unsaturated flow, 
standard practice in petroleum reservoir engineering is to introduce the following quantities: 

- liquid and gas formation volume factors, #/ = S/(p/), B g = B g (p g ); 

- solution gas/liquid phase Ratio R s = R s (/?/ ) ; 

- vapor water/gas phase Ratio R v = R v (p g ). 
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The explanation of formation volume factors and solution component/phase ratios, as used 
in the oil reservoir modelling, is as follows. Considering the volume AV[ es of liquid at 
reservoir conditions (reservoir temperature and pressure); when this volume of liquid is 
transported through the tubing to the surface it separates at standard conditions to a volume 
of liquid AVf d , and a volume of gas AV g ' d (coming out of the liquid, due to the pressure 
drop). At standard (i.e. stock tank) conditions, the liquid phase contains only oil component 
(here only water component) and the gas phase contains only gas component (here only 
hydrogen). Then applying conservation of mass, 

AVf es pi = AV^pf + AV s g ' d pf; 

and denoting the solution gas/liquid Ratio R s = AV g td /AVj" d , we may write: 

AV[ es Pl =AV? d {pf d +R s pf ), 
and the liquid phase mass density decomposition 

pf d +R s pf 

Pl= p ; (14) 



where B\ is the liquid formation volume factor, B\ = AV[ es / AV^ . Similarly, from AV g res p s 

V™/AV s g l 



and the gas formation volume factor B g = AV g es / AV g ' d , we get the gas phase mass density 



decomposition 



pf+R v pf d 
B 8 ; 



(15) 



where the vapor water/gas phase ratio R v = AVf d jAV g d . Now, in order to express the 
diffusion fluxes in terms of the new variables R S ,R V ,B/ and B g we have to rewrite the molar 
concentrations defined in 

i_ Wf cW _SiPf d ch JjPf_ cW _ s g R, P f d 

' M h B, ' ' M W B, ' * M h B g ' g M K B g 

These concentrations, defined in dl6l > may be slightly different from those previously defined 
in l|9]l if the gas at standard conditions is not composed only of hydrogen, and the liquid, also 
at standard conditions, is note made only of water. The phase molar concentrations, from 
d!6l >. are then given by: 



+c *--bA1¥ + -W ~ BZW il+FRv) > 



where F is given by: 



M h pf d 
M w p-l ld ' 



(17) 
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The component molar fractions in phase, are now defined as: 

xf = C -L 
ci 



(-'I 



Rspf 


R s 


R s p s g ' d + (M h /M w )pf d 


Rs + F' 


P\ td 


F 


pf ,d + {M W /M h )R s p s g ' d 


R s +F' 


p std 


1 


p s g ,d + (M h /M w )R v pf d 


l+FR v 


Rvpf 


FR V 



(18) 



* c g R v pf d + (M K /M h )p s g ' d 1+FR V ' 

Mass diffusive fluxes of components, defined in d 12b . depend on component molar fractions 
in phases; they have then to be rewritten. For instance, the mass diffusion flux of hydrogen 
in water, takes now the form: 

j*/pf = -cp^L^L^l(R s+ F)D l !VX" = -<p?L-!—D l tVR s . 
J//Ps p s g ' d BiM hy ' 1 1 B,R S +F' 

Other diffusion fluxes could be obtained by similar calculations, leading to the following 
formulas: 

<h h = t/ p f = -cp^L^—dlVR„ f!=jY/nf' d = 4>——^—DyVR s , (19) 
vi MiFg B,R S +F 1 " Vl J/m B,R S + F 1 " v ; 

S h =f/pf = <p^— - — D'lVR Vl A*=j»/o* , = -*5s. 1 D™VR V . (20) 

Vg *givg B g l+FR V g ' Vg Jj/P/ B g l+FR V s ' 

hr\h _ 



Remark 4 Let us recall, like in Remark [3j that the diffusion coefficients satisfy M h D 



Hi \ B\ j jj, g \ B g 

and the component fluxes, normalized by standard densities, are, 

<l> w = ^q; + ^q, + W + <i> h = j^di + ^ + + ^ 
Finally, we can write mass conservation Q, ((SJl in the following form: 

d ( Si RyS g \ , fA w\_ dzwi n sid 



M W D™, M h D g = M" D g ; and, like in CD, Li=h,wK = 1, for any phase a e {l,g}, in (m. 
The Darcy fluxes |Q} can now De rewritten in the following form: 

kn / pf+RsPf \ kr g ( pf+R v pf d 

q/ = -K— Vpi - — — g , q g = Vpg 



(21) 



*Si\B, + B) +diV {r) = ^ W/Pi ' (22) 

These equations have to be completed by equations $2$ and (fJJ. For instance, in i22[ and 
( 1231 ). we can take saturation and one of the pressures as independent variables; for example 
S g and pi, and for the other terms the following functional dependencies: 

Bi{pi),Bg(pg),R s (pi),R v (pg),Hi(pi),Hg(pg),kr w (Sg),krg(Sg). 

According to their definition, F, p g d , pf d andZ)' a (a 6 {l,g},i € {w,h}) are constants, and 
<P and K are depending only on space position. 
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2.3.1 Unsaturated flow 

Equations d22l > and d23l >. with saturation and pressure as unknowns, are valid if there is 
no missing phase (saturated flow); but if one of the phases is missing (unsaturated flow), 
equations and unknowns have to be adapted. There are two possible unsaturated cases, 
according to either the gas phase or the liquid phase disappears : 

1. Gas phase missing (Hydrogen totally dissolved in water): then we have S g = 0, Si = 1. 
Generalized gas phase Darcy's velocity is equal to zero since kr g (0) = 0. Independent 
variables are now pi, the liquid phase pressure, and R s , the solution gas/liquid phase 
Ratio; then we must write B; = Bi(pi,R s ) and ft = fij(pi,R s ), for < R s < R s (pi), 
where R s (pi) is equilibrium solution gas/liquid phase ratio. 

2. Liquid phase missing: then, S; = 0, S g = 1. Generalized liquid phase Darcy's velocity 
is zero since kr w (S g = 1) = 0; independent variables are now p g , the gas phase pressure, 
and R v the water vapor/gas phase ratio, become the new independent variables. However 
pressure p\ can be kept as independent variable since p g could be expressed through the 
capillary pressure law (3). We must also write B g = B g (p g ,R v ) and [i g = fi g (p g ,R v ), for 
< R v < R v (p g ), where R v (p g ) is the equilibrium water vapor/gas phase ratio. 

These above conditions can be summarized as 

S g >0, R s {p t )-R s >0, {R s (pi)-R s )S g = 0, (24) 
Si>0, R v { Pg )-R v >0, {R v ( Pg )-R v )S,=0. (25) 



Remark 5 In this last section, Phase Equilibrium Black-oil model, we did not take in account 
a possible interplay between dissolution and capillary pressure. □ 



2.4 Thermodynamical equilibrium Henry-Raoult model 

Another way of closing the system of equations 0, JQ, (O and $Q is to use phase thermo- 
dynamical properties for characterizing equilibrium . We use first ideal gas law and Dalton 
law, 

P g = Pg+P H g , (26) 

where p™ and p h are the vaporized water and hydrogen partial pressures in the gas phase; 
and 

p» = S^-RT, p\ = -\RT. (27) 

T is the temperature, R is universal gas constant and M w , M h are the water and hydrogen 
molar masses. 

Next, we apply Henry's and Raoult's laws which say that, at equilibrium, the vapor 
pressure of a substance varies linearly with its mole fraction in solution. In Henry's law the 
constant of proportionality is obtained by experiment and in Raoult's law the constant is the 
pressure of the component in its pure state. Here we will assume, for simplicity, that the 
quantity of dissolved hydrogen in the liquid is small; then these laws reduce to the linear 
Henry's law, which says that the amount of gas dissolved in a given volume of the liquid 
phase is directly proportional to the partial pressure of that same gas in the gas phase: 

pj'=H(T)M h p h g , (28) 
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where H(T) is the Henry's law constant, depending only on the temperature. For the liquid 
phase, we apply Raoult's law which says that the water vapor pressure is equal to the vapor 
pressure of the pure solvent, at given temperature, multiplied by the mole fraction of the 
solvent. The water vapor partial pressure of the pure solvent depends only on the temperature 
and therefore is a constant, denoted here by p"(T), so we have from definition Q 



p»=p«(T)Xr=p;{T)—£-^. (29) 



_p7 

' pf + (M W /M h )p[ J 

Further on, we can include in formula \2% the presence of capillary pressure, by using 
Kelvin 's equation (see [4|) which gives: 

w = -r/jn P± e -M-*'p c /(RT Pl ) (3Q) 

If now, to equations l|26ll-(|28ll and d30b we add the relations 

pf + pf = p h p h g +p» = p g , (31) 
and the water compressibility, defined by 

n std 

p ^4o ; <32) 

then, we have 8 equations: l !26l l. d27lh 9. I l28l l. d30b . Blti 9 and ( 132b and 10 unknowns: 

Pi , Pg , Pg , Pg , Pi , pf , Pi , Pg , Pg , Pg ■ 

We may, for instance, parametrize all these 10 unknowns, by the two phase pressures pi 
and p g . For this we should combine the Henry law l !28l >, the Raoult- Kelvin law < |30b and ( |26t 
leading to the system of two equations for p h g and p g : 

„ w = ff v ( T) — e -M" (p s -p,)/(RT( P }«+M f 'H(T)p h g ) 

h Pgy > p* + (MwH(T))pl 

Pg=Pg+Pg- 

It is easy to show that this above system of two equations has a unique solution pg,pg > 
for any p g > p g '(T); and solving this system we obtain 

Pg=f(Pg,Pl), Pg = g{Pg,Pi)- 

By d27|li 2 we can then write p g , p g and p f as functions of p/ and and by d28b and d32| i 
we can finally express p[\ pf and p\ as functions of the phase pressures. 

Remark 6 The gas phase will appear only if there is a sufficient quantity of dissolved hydro- 
gen in the liquid phase, and this quantity is exactly the dissolved gas quantity, at equilibrium, 
given by Henry's law. But when the dissolved gas (hydrogen) quantity is smaller than the 
quantity of hydrogen at equilibrium, then the Henry law does not apply; and S g is then equal 
to zero and could not be taken as unknown. In this situation, instead of saturation, we may 
take pj 1 as independent variable. We notice that Henry-Raoult model based on thermody- 
namical equilibrium leads to a similar concepts as the ones developed for establishing Black 
Oil model for reservoir modeling. □ 
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Remark 7 In the Henry-Raoult model, if there is no vaporized water, p g = pf, the criteria 
for non saturated flow is simple and reads 

5^<ft = w+*(°). (33) 

if there is a threshold pressure, i.e. p c (0) ^ 0. The same criterion can be used also if there 
is vaporized water, as long as the pressure in porous medium is much larger than vaporized 
water pressure. 

Remark 8 In Henry-Raoult model, above, we were assuming for simplicity, no complete 
evaporation of the water. 



2.4.1 Comparison of Equilibrium models 

We will now consider in more details a special case, where we may, from the Henry-Raoult 
model, get back the Black oil model. In both models we will then use diffusive fluxes devel- 
oped in Section [23] For simplicity, we will neglect in ( |29b influence of capillary pressure 
and solved gas on water vapor partial pressure, leading to p™ = p g (T) being a constant. 
Then we have from (HU, (1321 and (f3TT> 



,/, „„■ _ (f ,-,,„, V , _Pf_ = Pf d + B,(p,)H(T)M h (p g -p^T)) 
Vg B,(Pi) MPD 



p, = pf + pf = H(T)M»p'l + = 2 D _/„_, v " ^ , (34) 

and also from ( 127b . 



p g = p h g +p- = ^( Pg - p»(T)) + ^%(T). (35 ) 



Therefore, equations i34l and d35b can be written in a form similar to d 14b-d 15b by defining 
the gas formation volume factor B g , solution gas/liquid phase ratio R s and vapor water/gas 
phase ratio R Y , from above thermodynamic al relations: 

R s =B l {p l ) t ^^{ Pg -p^{T)) (36) 
"g 

RTpf _ 1 p»(T) 

Bg ~ M"(p g -p-(T)Y Rv ~Fp g -p-(T)' (3?) 

where F is given by ( 1171 ). 

Compared to the Black-Oil model, Section [231 here, it is clear from d36l l that solution 
gas/liquid phase ratio R s depends on p\ and p g and not only on p\. 

Remark 9 For the case without any water vapor, the corresponding thermodynamical model 
is obtained simply by taking p\ v {T) = in d36l l. ( 1371 ) and R v = in equations ( 122b and l !23b . 
In the same way, if we neglect dissolved hydrogen, the corresponding thermodynamical 
model is obtained simply by taking H(T) = in ( 136b . ( 137b . and R s = in ( 122b and l !23b . □ 
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2.4.2 Model assuming water incompressibility and no water vaporization 

In the Henry-Raoult model we assume now that the water is incompressible and the water 
vapor quantity is neglectful; i.e. the gas phase contains only hydrogen, then (see Remark[7} 
p h g = p g and p H g = pJ{T) = in 03 leading to R v = 0. Formulas OSK 03». could be 
rewritten as 

B t = 1, R s = C h p g , — = C v p g ; (38) 

Br, 



where we have denoted 



_ H(T)M h M h 
Qx ~ pf ' Cv ~Rfp~f ; (39) 

and equations d 1 4b— <l 15b become: 

p,(R s ) = pf d +R s pf, Pg(pg) = Q.pfpg. (40) 

Similarly to constant F defined by d 1 Vb , we use the density ratio 

std 

G=^7. (41) 



and equations d22l >. d23b reduce to 



*^+div(q, -I ) - ■/ "■/>;"'. (42) 



q/ 



* ^ (5A + C vPg S g ) + div (j^q, + C vPg qg + j) =& h /pf, (43) 

Vp, - (pf d +R s pf)g) , q g = -K^ (v Pg - C„pf Pg g) , (44) 

<f>S/,F 



D ? V^. S , (45) 



where we have denoted 0* = J and used Remark[3]to get FD* = GZ)J'\ from which it follows 
W = -J/G. 

Here there is only hydrogen in phase gas, pi = p g and Henry's law assumes thermody- 
namical equilibrium in which quantity of dissolved hydrogen is proportional to gas phase 
pressure. In saturated case (where two phases are present) Henry's law reads R s = C/,/? g , and 
we can then work with variables p g and Si in equations (I42t-l|45ll. 

When the gas phase disappears, the gas pressure drops to the liquid pressure augmented 
by entry pressure, p g = pi+ p c {0), the liquid can contain any quantity of dissolved hydrogen 
p'g 1 between zero and C h p s g td (pi + p c (0)) = H(T)M h (pi + p c (0)), from Henry's law (see <[28j 
and definition d39l>). 

And then, when the gas phase is absent (one of the unsaturated cases) S g = 0, we will 
replace saturation, as we did in the Black-Oil model for the same unsaturated case in Sec- 
tion |2.3.l1 by a new variable R s (see definitions l !28l >, d36l l, l!39l>), such that 

is the mass density of dissolved hydrogen in the liquid phase. 
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Assuming at standard conditions the gas phase contains only hydrogen (hydrogen com- 
ponent masswgas phase mass) and the liquid phase contains only water, with water in- 
compressibility and no vaporized water, we see that, in definition I l36t . R s = p/'/Pf"' — 
AVg' d I AVf d which is exactly the definition given in the Black-Oil model in Section [2~3l 

The physical meaning of R s then stays the same either in Henry-Raoult or in Black-Oil 
models even in the unsaturated case. Moreover, from mass conservation law it follows that 
the dissolved hydrogen mass density R s p*' d is continuous when the gas phase vanishes and 
this can be expressed as an unilateral condition: 

0<Sg<l, 0<R s <C h p g: S g (C h p g -R s ) = 0. 

In unsaturated region, where S/ = 1 (that is S g = 0) we replace variable Si by R s , and 
equations d42b — d45 b degenerate to: 

div(q,-ij) = & w /pf d ; (46) 

+div(fl s q, + j) = & h lpf; (47) 

q, = -KA,(1) (Vpi - ( P f d +R sP f)g) ; (48) 

J=--—D'IVR S . (49) 

K. + 1 



2.5 Saturated/unsaturated state, general formulation 

Finally in the saturated region we used p\ and S g as variables in d42b-(|45ll but in unsaturated 
region we should use other variables, p\ and R s , in d46b — d49 b - In order to avoid the change of 
variables and equations in different regions as above we prefer to introduce a new variable 

X = (l-S g )Rs + C v p g S g ; (50) 

in view to make equation d43b parabolic in X. 

This new variable X is well defined both in saturated and unsaturated regions. Moreover, 
from <[3Sj and ([28]l, R s = pf/pg d ; from <|40j C v p g = Pg/p s g ' d and since p g = this new 
variable X is a "normalized total hydrogen mass density": X = (5/p/ 1 + S gP g)/ Pg " d . It 
is easy then to see that parabolicity is possible only if we take liquid pressure p\ as other 
independent variable. Knowing that in saturated case, S s > 0, R s = Ckp g from Henry's law, 
we may write X defined in d50b as: 

x= {(C h (l-S g )+C v Sg)( Pl+Pc (S g )) ifS g >0 
[R s ifSg = 0. 

Since for capillary pressure defined as a function of gas saturation we have p' c {S g ) > 0, and 
since usually, like for hydrogen, in ( |39b CO = C v /Ci, > 1 we get the following bounds: 



a(Sg)=C h (l-S g )+C v S g e[C h ,C v }, a'{S g )=C v -C h =C A >0. (52) 
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Then for any S g > from i5H 
dX 



ds g 



= C A ( Pl +p c (S g ))+a(S g )p' c (S g )>0, 



and therefore for each p\ > we can find inverse function S g = S g (pi,X) which satisfies 

d S 

> for S g > 0. 
After taking derivatives with respect to pi and X of d51| l we obtain: 

dS g _ a(S g ) 2 x(p,,X) dS g _ a(S g )x(pi,X) 

dp, C A X + a{S g ) 2 p' c {S g Y dX C A X + a{S g ) 2 p' c {S g y 

where X ( P i,X) ls characteristic function of the set {X > Ci,(pi + Pc (0))}. In j53K we remark 
from l l52t and definition ® that dS g /dpi < 0. 
Let us note that property 

P ' c (S g = 0) = +00, (54) 

of van Genuchten p c functions, leads to continuity of the two above partial derivatives since 
we have 

dS s ,. dS e 
hm -r-^ = lim ^ = 0. 

s g ^o dpi Sg^O dX 

We now introduce auxiliary function N(p/,X) defined as 



C A X 

C A X + a(S g ) 2 p' c (S g y 



N (P'> X ) = T> v , „,c X2 77* a (Pi,X) G [0, 1) 



which verifies 



X in,xn P :<s g )^.=N{p h x), p ' c{Sg )^ = L-^l x ( Pi .xi ,55. 



Note that function N( P i,X) is continuous under condition 

Darcy's fluxes in d44b and diffusive flux d45b now take the form, with Pg = P i + Pc (S), 

q, = -KX,(S g ) (V PI - (pf d + R s ( Pl ,X)pf)gj (56) 
q g = -KX g (S g ), (Vpi+V Pc (S g )-C v pf Pg ( Pl ,X)^ , (57) 
,D';VR s ( Ph X); (58) 



R S ( PU X)+F 

where S g is a function of P \ and X and where R s , defined as in d3 8 b . can now be expressed 
in both saturated and unsaturated region as a function of new variables P i and X: 

R s (pi,X) =mm(C hPg ( P i,X),X), p g (pi,X) =pi+p c (S g (pi,X)). 
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After expanding the capillary pressure gradient in d57l l and the gradient of R s in d58b . we 
may write: 

q g =-KX g (S g ) (ll+p^^Vn+p'^-^VX-Cvpfpsip^g 



<P(\-S g (p,,X))F ^ dS g 
J = R s ( Ph X)+F VlW+Pc^^Pl 

0(1 -S g (PuX))F dS, *F_ 

R s (p h X)+F ' huys 'dX X + F n AJ 

From d2 lb water component flux <j) w and hydrogen component flux <p h , with the assumptions 
of incompressible water and absence of water vapor, reduce to 

f r = q/ + = -(A u V Pl +A L2 VX + S,), (59) 
$ h =R s qi+C vPg q g + tf = -(A 2A Vpi+A 2 , 2 VX + B 2 ), (60) 



where the coefficients Aij(pi,X) and Bi(pi,X) are given by the following formulas (note 
b'l = J and 



that from © $ = J and (j)) 1 ' = -J/G): 



A n ( Pl ,X) =KX,(S g ) - ^—^-D';c h N, (61) 
- , . <P(l-So)F 1 -N h 

A\ lipuX) = - — ^ —D,C h , (62) 

' u ; {R S +F)G a(S g ) 1 ' v ; 

A 2 .i (p/,X) =KA,(^)^ + KA g (^)C,. Pg jV+ ^ ~J g J F D';c h N, (63) 

A 2 , {m ,X) =^S^C vPg + ^^^C h , (64) 

=-KA / (5 g )[p/' rf +,R J pf ]g, (65) 

2fe(p,,X) =-KX,(S g )R s [ P f d +R sP f]g-KX g (S g )Cf, P f P 2 g g-, (66) 



where we have used d55l l and A. g (S g )(l — %) =0. 
Equations J42b-d45|l become: 

-4>|^ -div (A M V W +A I , a VX+B 1 ) - = (67) 

4>^-div(A2.iV w +A 2 ,2VX + 5 2 ) =& h /pf. (68) 



The gain in this form is that d68l > is a parabolic equation for X since 

A 2 . 2{PI ,X K ■ 4 =K^ ■ |A,(S g ) i^C^ + i^DfCH |» 

is strictly positive in the whole domain if the diffusion and capillary pressure are not ne- 
glected. 

If we eliminate diffusive terms from equation d67| l (the "pressure equation") by forming 
an equation for total flow <j> tot , defined in l l69l l. that is summing equation d68l and equation 
| |67| |. we obtain 

<l> tot = G4> K + (t> h = (G + R s )q,+C vPg q g = -(A M V W +A 1|2 VX+B 1 ), (69) 
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where G is the density ratio defined in d4 1 b . and coefficients A\ \ ,A\ 2 and B\ are given by: 
A 1A (p,,X) =GA\j +A 2J = KXi(S g )(G + R s )+Kk g (S g )C v p g N, (70) 
A hZ ( Pl ,X) =GA 12 +A 2 , 2 = KX g (S g )^-C vPg) (71) 

Bi(j>i,X) =GBi +B 2 = -Kh{S g )(G+R s )[pf d +R s pf}g-KX g (S g )Ctpfp z g g. (72) 

Now the "pressure equation" (67l is transformed to: 

dS d v dS dX 

-G<Pj^--^--div{A lA V Pl +A 1 , 2 VX + B l ) + *(1 - G-=£)-£ (73) 

= G& w /pf d + & h /pf. 

With this last formulation, we see that the "pressure equation" l !73l > is parabolic/elliptic equa- 
tion in pi since 

Au(pi,X)£ = K| ■ |XK^)(G+R,)+K| ■ !A g (S g )C v /VV, 

is strictly positive, independently of presence of diffusion or capillary forces, and the coef- 
ficient in front of dpi/dt is positive since, as we remarked at d53l >. dS g /dpi < 0. 

Finally, the transport of water and hydrogen is described by differential equations fl73t 
and d68l > which are rewritten here, using ( 160b and d69l >. in the form 

<Pj t (X-GS g ( Pl ,X))+div($ lol ) = GJ?™/pf d + & h /pf, (74) 

G^+div^^^/pf, (75) 

where the fluxes are given by d60l > and l |69b , while the coefficients are given by formulas 
(f61~b— J66t and 

2.5.1 Boundary conditions 

Equations d74b and d75b . given in porous domain Q , must be complemented by initial and 
boundary conditions. Following ([3 D, we assume that the boundary dQ is divided in several 
disjoint parts: impervious, inflow and outflow boundaries. We present now a set of standard 
boundary conditions on each of these boundary parts. 

- On impervious boundary F mp we take Neumann conditions: 

(j, tot . v = and <j> h -v = 0. (76) 

- On inflow boundary when pure water is injected, we impose for hydrogen component 
X = and for liquid pressure either p\ = p\ ,„ or <j) lot ■v = Q l j. 

- On inflow boundary, when pure gas is injected we have (j) H ' ■ v = and we can impose 
for the pressure, total injection rate which is then equal to gas injection rate: 

* t0( -V = G5 = **.V. (77) 

- On the outflow boundary, when liquid is displaced by gas we have possibility to impose 
for gas phase X = and for liquid pressure p\ = pi. out , only before gas reaches this 
outflow boundary (breakthrough time). Alternatively, with the same Dirichlet condition 
for liquid pressure, for gas saturation we can set either Neumann condition VS g ■ v = 0, 
or Dirichlet condition S„ = 0. 
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2.6 Model without capillary pressure, diffusion or gravity 

When capillary pressure is neglected, we have only one pressure p g = p\ = p, then definition 
d51l > of variable X simplifies to: 

x= Uc h {\-s g )+c v s g ) P ifS g >0 

[R s HS g = 0, 



<P^(X-GS g ( P ,X)) + dtv((t> lol ) = GJ? w /pf d + & h / P s J d (79) 



and partial derivatives d53 b can be calculated explicitely , 

^ = —pr-iXiP,*), |S = J-X(P,X): 
op Cap aX C^p 

where %{p,X) is the characteristic function of saturated region, as before. 

When we neglect capillary pressure, and also diffusive and gravity fluxes, system J74b , 
d75l l reduces to the following two equations 

dt" 

4>^ + div(/ A (p,X)tfy) =J?"/pf , (80) 

where 

$ m =-A of (p,X)KVp, $ h = f(p,X)4> tot . (81) 
Total mobility A tot (p,X) and hydrogen fractional flow functons f h (p,X) are defined by: 

A tot (p,X)= h (S g ) (G + R S ( P ,X))+ A s {S g )C vP , 

A h (p,X) = ?n(S g )R s (p,X) +A g (5,)C vA f h ( P ,X) = 

System (I79M8 1 ll is very close to immiscible two-phase system (see for instance |3 1); the 
only difference in hydrogen transport equation d301 ) is the presence of R s (p,X) in fractional 
flow function f h (p,X), and if we write pressure equation l !79t in the form 

G ^^2 If + div (^) + *(1 - §^)^f = GJ? w /pf d + J? h /pf, (82) 

then we see that the presence of dissolved gas introduces additional "source" term in total 
flow equation l[79]>, namely <5(1 - G%l{C A p))dXjdt. 

Remark 10 Considering the boundary conditions defined in the previous section; on the 
boundary part where there is hydrogen outflow, to impose a Dirichlet condition on X will 
lead to a boundary layer. 



3 Numerical simulations 

This section presents two test cases and their simulations using the new formulation given 
by system d V4b . d75l l. Both test cases are not build to correspond with particular real sit- 
uations but rather to illustrate the gas appearance phenomenon. Assuming horizontal two 
dimensional problems gravity effects are neglected in both cases. The first test case is a one 
dimensional like situation where hydrogen is injected through an inflow boundary and the 
second test case is a two dimensional situation where hydrogen is injected via a volume 
source term. 
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3.1 Setting test cases 
3.1.1 Physical data 

In the two test cases we consider the same isotropic porous medium with a uniform absolute 
permeability tensor K = k where k is scalar and a uniform porosity <P. The capillary pres- 
sure function, p c , is given by the van Genuchten model (see 1 12]) and relative permeability 
functions, kr\ and kr g , are given by the van Genuchten-Mualem model (see [ 12| and 1 13 1). 
According to these models we have : 

Pc = Pr (S u 1/m - l) ^ , kn = VWe (l - (1 - S)i m )) 2 and kr g = s/T=ST e (l - s]'" 1 ) *" 
with Si e = and m = 1 

I Sir- - Sgr tl 

where parameters P r , n, Si r and S gr depend on the porous medium. Values of parameters 
describing the considered porous medium and fluid characteristics are given in Table [T] 
Fluid temperature is fixed to T = 303K. 



Porous medium parameters 


Fluid characteristics 


Parameter 


Value 


Parameter 


Value 


k 


5 W- m 


m 1 


D>1 


3 10- y 


m l js 


<P 


0.15 


(-) 


W 


1 10~ 3 


Pa.s 


Pr 


2 10 6 


Pa 




9 10- 6 


Pa.s 


n 


1.49 


R 


H(T = 303K) 


7.65 10" 6 mol/Pa/m 3 


Sir 


0.4 


(-) 


M, 


10- 2 


kg/mol 


Sg r 





(-) 


M g 


2 10~ 3 


kg/mol 








Pf d 


10 3 


kgjrr? 








Pf 


8 10~ 2 


kgjrr? 



Table 1 Values of porous medium parameters and fluid characteristics 



3.1.2 Test case 1 

In the first test case we consider the domain Q 1 = [0m ; 200m] x [—10m ; 10m] with an 
impervious boundary F^ mp = [0m ; 200m] x {—10m, 10m}, an inflow boundary F^ n = {0m} x 
[—10m ; 10m] and an outflow boundary F} ut = {200m} x [—10m ; 10m]. The following 
boundary conditions are imposed : 

- <j> tot ■ v = ty h ■ v = on the impervious boundary r^ np , 

- /of • v = <p h ■ v = on the inflow boundary , 

- X = and pi = pj out on the outflow boundary F} ul 

where Q h d and p\ mt are constant scalars. Source terms are fixed to zero = ^ w = 0). 
Initial conditions are X(t = 0) = and p/(t = 0) = p\ out on Q} . The boundary parameters 
are fixed to Q h d = 1.5 10~ 5 m/years and p\ mt = 10 6 Pa. 
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3.1.3 Test case 2 

In the second test case we consider the domain Q. 2 = [Om ; 200m] x [— 100m ; 100m] with an 
outflow boundary T 2 U = dQ 2 and where B\ = [90m ; 1 10m] x [—10m ; 10m] is the support 
of hydrogen source term. On the outflow boundary F 2 ul we impose X = and p\ = pf out . 
Source terms are defined by ,jfi h = F?x B 2 an d ^ w = 0- Initial conditions are X(t = 0) = 

h 

and pi(t = 0) = pf oM on Q 2 . Here pf oM and are constant scalars fixed to p 2 oul = 10 6 Pa 
and Fl = 8 10~ 13 kg/m 3 /s « 2.5 10~ 5 kg /m 3 /year. 

3.2 Numerical results 

System ( I74l l. d75l > is a coupled nonlinear partial differential equation system. Numerical sim- 
ulations use an implicit scheme for time discretization, a finite volume scheme (using Multi 
Point Flux Approximation) for space discretization and a Newton-Raphson like method to 
solve nonlinearities. All computations are performed with the Cast3m software (see 1151 ). 

In both cases we present, at several times, spatial evolutions of the liquid pressure, the 
total hydrogen molar density and the gas saturation along the line ££ cut = [0m ; 200m] x 
{0m}. Computations are performed since the time t = up to the stationary state. 

3.2.1 Results and comments 

Results of test case 1 are plotted on figures [TJ [2] and [3] The total hydrogen molar density 
-gj-X (FigfTJ, the liquid pressure pi (Fig[2]> and the gas saturation S g (Fig[3jl are plotted 
at times t = 1 10 4 , 2.5 10 4 , 5 10 4 , 1.1 10 5 , 2.5 10 5 and 5 10 5 years. Results of test case 2 

are plotted on figures |4l [5] and [6] The total hydrogen molar density -nwX (Fig|4}, the liquid 
pressure pi (FigO and the gas saturation S g (Fig|6]l are plotted at times t = 50.1, 125, 355, 
2820, 2 10 4 and 10 5 years. 

In both cases, we can identify three characteristic times : at t = T\ the gas phase ap- 
pears; at t = T2 the maximum liquid pressure is reached; at t = Tt, the system is close to the 
stationary state. For test case 1, we have 

T\ « 2 10 4 years, T2 ~ 1.1 10 5 years and T3 w 5 10 5 years; 

for test case 2, we have 

T\ w 90 years, T2 « 355 years and T3 w 10 5 years. 

Global behaviors of both cases are similar and can be summarized as follow : 

- For < t < T\ : only total hydrogen density increases while liquid pressure and gas 
saturation stay constant; during this stage X < C/,pi and all the domain is saturated in 
water (S g = 0). 

- From t = T\,X> Chpi in a part of the domain meaninig that gas phase exists (S g > 0) 
in this part. 

- For 7] < t < T2 : while gas phase appears, liquid pressure increases and a non zero pres- 
sure gradient appears what corresponds to a fluid displacement according to the Darcy- 
Muskat law. Total hydrogen density and gas saturation increase and the unsaturated area 
grows. 
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- For ?2 < t : while total hydrogen density and gas saturation continue to increase, liquid 
pressure and pressure gradient decrease. When t — > °°, the system reach a stationary state 
where saturated and unsaturated areas coexist and liquid pressure gradient is null. 



16 



-- 14 
>> 

w 

£ 12 < 

ro 10 



* 5 



i. < 



A 



80 100 120 

abcissa (m) 




W 1 W ■ 



t=1 .00E+04 
t=2.50E+04 
t=5.00E+04 
t=1.10E+05 
t=2.50E+05 
t=5.00E+05 



A 



140 160 180 



&4 



Fig. 1 Test case 1 : spatial evolution along the line _S? C „, of the total hydrogen molar density -j^-X at several 
times t (in years) 
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Fig. 2 Test case 1 : spatial evolution along the line ££ cut of the liquid pressure p\ at several times t (in years) 
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Fig. 3 Test case 1 : spatial evolution along the line _S? C „, of the gas saturation S g at several times t (in years) 
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Fig. 4 Test case 2 ! spatial evolution along the line J£ cu 
times t (in years) 



of the total hydrogen molar density -j^X at several 



4 Concluding remarks 

From balance equations, constitutive relations and equations of state, assuming thermody- 
namical equilibrium, we have derived a model for describing underground gas migration in 
water saturated or unsaturated porous media, including diffusion of components in phases 
and capillary effects. In the last part of this paper, numerical simulations on simplified sit- 
uations inspired by the "Couplex-gas" benchmark |14], show evidence of its ability : - to 
describe gas (hydrogen) generation and migration - and to treat the difficult problem, as it 
appeared in the results of "Couplex-gas" 1 14|, of correctly simulating evolution of the un- 
saturated region, in a deep geological repository, created by gas generation. A forthcoming 
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Fig. 5 Test case 2 : spatial evolution along the line l£ cut of the liquid pressure p\ at several times r (in years) 
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Fig. 6 Test case 2 : spatial evolution along the line Jz? c „, of the gas saturation S g at several times t (in years) 



paper will be devoted to the use of this model for solving the "Couplex-gas" benchmark 1 14 1 
and other 3-D situations of gas migration in water saturated or unsaturated porous media, 
including the design a of numerical test cases synthesizing the main challenges appearing in 
gas generation and migration. 
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